Monte Carlo study of a two-dimensional quantum ferromagnet 



On 
On 
Os 



c3 

CZ2 



c3 



i 

C 

o 
o 



> 

00 

m 
o 

On 
Os 

-i— > 
03 



i 

-a 
c 

o 
o 



Patrik Henelius a , Anders W. Sandvik 6 , Caxsten Timm c and S. M. Girvin 
a National High Magnetic Field Laboratory, Tallahassee, FL 32310 
Department of Physics, University of Illinois at Urbana- Champaign, Urbana, Illinois 61801 
: Institut fur Theoretische Physik, Freie Universitat Berlin, Arnimallee 14, D-14195 Berlin, Germany 
Department of Physics, Indiana University, Bloomington, Indiana 4^405 
(February 1, 2008) 

We present quantum Monte Carlo results for the field and temperature dependence of the magneti- 
zation and the spin- lattice relaxation rate 1/Ti of a two-dimensional S — 1/2 quantum Heisenberg 
ferromagnet. The Monte Carlo method, which yields results free of systematic errors, is described 
in detail. The high accuracy of the calculated magnetization allows for stringent tests of recent 
approximate analytical calculations. We also compare our results with recent experimental data for 
a v = 1 quantum Hall ferromagnet, which is expected to be well described by the Heisenberg model. 
The dynamic response function needed to extract 1/Ti is obtained using maximum-entropy analytic 
continuation of the corresponding imaginary-time dependent correlation function. We discuss the 
reliability of this approach. 



I. INTRODUCTION 

The magnetization of a two-dimensional electron gasi 
the quantum Hall regime has recently been measured a 
The realization that this system represents a novel itiner- 
ant_fcrromagnet has motivated several theoretical stud- 
ies.fTTl The system at filling factor v = 1 should be 
well described by a two-dimensional Heisenberg ferro- 
magnetu, and in a recent publication we presented analyt- 
ical Schwingerj-boson and numerical Monte Carlo results 
for this modelE The ferromagnetic Heisenberg model has 
been much less studied than the antiferromagnetic model, 
probably because the ground state and lowest excitations 
are known exactly. At finite temperatures there is, how- 
ever, no exact analytic solution, and in our first paper we 
gave a detailed discussion of the calculation of the lead- 
ing 1/N correction to the Schwinger boson mean-field 
theory. In this second paper we want to give a detailed 
description of the Monte Carlo calculation and present 
results for the NMR relaxation rate 1/Ti. The Monte 
Carlo technique used is highly efficient and suffers from 
none of the common systematic errors. The magnetiza- 
tion results are accurate to a relative statistical error of 
about 10~ 4 and to this accuracy they are also free of 
finite size corrections. This allows for stringent tests of 
the analytical results. In order to estimate 1/Ti we calcu- 
late spin correlation functions in imaginary time and con- 
tinue these to real frequency using the maximum entropy 
method. We discuss how the results of this approach de- 
pend on whether real- or momentum space correlation 
functions are analytically continued. 

We want to emphasize an important technical detail 
that makes the sampling particularly efficient: the ex- 
ternal field is chosen in the x-direction, perpendicular to 
the z-direction, in which the basis states are expressed. 
This automatically causes the simulation to become free 
of systematic errors, even though only local updates are 
used. Furthermore, it enables easy access to observables 
involving both diagonal and off-diagonal operators. 

In Sec. II the stochastic series expansion Monte Carlo 



technique is briefly reviewed and applied to the two- 
dimensional ferromagnet. The maximum entropy tech- 
nique used to numerically perform the analytic continu- 
ation to real frequency is described in Sec. III. Results 
for the field and temperature dependence of the magne- 
tization are given in Sec IV. A, and the NMR relaxation 
rate results are presented in Sec IV. B. Our results and 
conclusions are summarized in Sec. V. 



II. QUANTUM MONTE CARLO 
A. Introduction 

Over the last few decades, several methods for quan- 
tum Monte Carlo calculations have been developed. The 
most common methods stochastically sample the config-i 
urations of the world-line path integral of the system.EI 
These methods traditionally suffer from several system- 
atic errors.113 The Trotter discretization of the imaginary 
time introduces a systematic error that decreases with 
decreasing "time slice" width At. In principle exact re- 
sults are obtained by— scaling to At — » 0. However, with 
standard techniques, til the efficiency, of the simulation de- 
creases rapidly as At is decreasedJ13 and it may therefore 
be difficult to obtain results completely free of systematic 
errors. Each configuration can furthermore be labeled 
by a topological "winding number" . The acceptance rate 
for moving from one winding number sector to another 
gets extremely small as the system siae increases, thereby 
making the simulation non-ergodiclia In addition, at low 
temperatures it can be very hard to change the number 
of particles in fermion or boson simulations, or the total 
magnetization for spin systems, hence restricting simula- 
tions to a canonical ensemble. 

Recently there has been much progress in resolving 
these technical problems, making it much easier to obtain 
results that are exact within statistical error bars. Non- 
local "loop cluster" moves (analogous to cluster updates 
for classical systemsE3) have been developed for several 
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modelsBO and can significantly reduce the autocorre- 
lation times of the simulations. They enable sampling 
of all winding number sectors and all magnetizations. In 
addition these moves have been formulated in continuous 
imaginary timc,ll3 hence eliminating the Trotter error and 
making the simulation completely free of systematic er- 
rors. A related method was recently formulated starting 
from the standard perturbation expansion in the interac- 
tion representation.^ Within this formulation, new local 
moves that share some of the advantages of the loop al- 
gorithms was introduced. 

A different approach to quantum Monte Carlo simu- 
lation isr-the, so called stochastic series expansion (SSE), 
method liZHia (a generalization of Handscomb's methodt£l 
to a much larger class of models) which does not use the 
Trotter decomposition as a starting point, but instead 
samples matrix elements of the exact Taylor expansion of 
the density matrix e~^ H . The Trotter error is thus auto- 
matically avoided. There are strong formal relationships 
between the SSE formulation and the continuous-tip/: 
path integral, which have been discussed elsewhereEl 
We will here demonstrate that when applying the SSE 
method to a ferromagnet in a magnetic field, purely lo- 
cal moves are sufficient to make the calculation ergodic 
in the grand canonical ensemble. 

Below we will first give a short general overview of the 
SEE method, and thereafter give a detailed description 
of the application to the Heisenberg ferromagnet in a 
magnetic field. 



B. Stochastic Series Expansion 

The thermodynamic expectation value for an operator 
A, at inverse temperature /?, is 



(A) = \ = -Tr{Ae-P H \. 



(1) 



Assume that the Hamiltonian consists of M terms that 
need not commute, and may be of any order in field op- 
erators: 



(2) 



A minus sign has been pulled out of the sum for conve- 
nience. The density matrix exp(— (3H) in the expectation 
value Eq. (pi) is Taylor-expanded, 
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(3) 



where S n denotes an index sequence hi, . . . , k n with 
1 < h < M, and Yli =1 H ki is an operator string of length 
n. If the above trace can be analytically calculated the 
expectation value can be calculated by importance sam- 
pling in the space of index sequences. In the original 



Handscomb's methodB a spin-1/2 system is considered, 
for which the trace can bp-evaluated analytically. In the 
more general SSE methodic a complete set of states {\a)} 
is introduced to calculate the trace, and hence a larger 
class of problems can be treated: 
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(4) 



a n=0 S n 



i=l 



In order to calculate the expectation value of the operator 
A we assume that a function A(a, S n ) exists such that 
Eq. (0) can be re-written as: 



(A)= 



E a En= Es n A ^Sn)W(a,S n ) 

£a£~o£ Sn W(a,S B ) 
(A(a,S n )) w , 



where the weight function W(a, S n ) is given by 



a 

W(a,S n ) = tL.{ a \Y[H ki \a). 



(5) 



(6) 



We will assume that W(a, S n ) is positive definite, which 
normally is the condition for a stochastic evaluation of 
(||) to be feasible. One can then carry out a random 
walk satisfying the detailed-balance principle in the space 
{\a) ® {S n , n = 0, . . . , oo}}, using W as a relative prob- 
ability. One important condition for such a procedure to 
be feasible in practice is that the operators Hi in Eq. (||) 
must be "non-branching", i.e., the application of any Hi 
to a single basis vector should yield a single basis vector 
(not a linear combination of them), 



Hi\a) = h(i, a)\a'), 



(7) 



so that the weight factor (||) can be easily calculated. 
A scheme can then be constructed in which first an 
arbitrary allowed operator string and state are chosen. 
Thereafter relatively simple updating steps ( "moves" ) are 
performed that change the number of operators (n) in the 
string, the individual operators within the string (thereby 
changing S n ), and the state The acceptance prob- 
abilities for the moves are chosen so that the detailed- 
balance principle is satisfied, using, e.g., the standard 
Metropolis or heat-bath algorithm. 

The Taylor expansion may appear to be a high- 
temperature expansion, but in principle terms up to 
n = oo are included, and the expansion is equally valid 
at any temperature. For a finite-size system at a finite 
temperature only terms of finite length contribute signif- 
icantly to the trace, and importance sampling includes 
all relevant terms. This can be compared to the Taylor 
expansion of the exponential of a simple number, where 
the factorial n! in the denominator eventually wins over 
the numerator. An actual distribution for the order of 
the series in a typical simulation appears to be close to 
a normal distribution; see Fig. [I]. As will be derived be- 
low, the average length of the operator string, (n), equals 
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FIG. 1. Distribution of the order of the Taylor expansion 
for a 4 x 4 Heisenberg ferromagnet at temperature T/J = 2.0 
and magnetic field h/J = 1.0. 



f3E, where (3 is the inverse temperature and E is the to- 
tal energy. At low temperatures we thus see that the 
average operator-string length is inversely proportional 
to the temperature and proportional to the number of 
sites. The computational time required for one MC step 
is proportional to the inverse temperature times the sys- 
tem size. The computational cost to achieve a given ac- 
curacy is, however, often offset by the fact that as the 
length of the operator-string is increased it also contains 
more information about the system. This will become 
clear when we discuss how to measure various expecta- 
tion values below. 

Next, we need to find the function A(a, S n ) for differ- 
ent cases of interest. First we look at some important 
features of the states and introduce some notation. It is 
convenient to introduce states \ct(p)) that are obtained 
by propagating |a) with p (p = 1, . . . , n) of the operators 
in Eq. (p|); 



Y{p))=ri[H ki \a), 



(8) 



where r is a normalization factor. With this definition 
|a(0)) = \a). Using this notation the matrix element in 
the weight function is 



(9) 



UtiHh, a(i-l)], if \a(n)) = |a(0)) 
0, otherwise. 



If the index sequence S n is cyclically permuted p times we 
obtain /c p +i, . . . , k n , fei, . . . , k p , which is denoted S n (p). 
Since W(a, S n ) = W(a(p), S n (p)) it follows that 



{A) = ( 7 ^-iTAla(p),S n (p)} 

p=0 



w 



(10) 



In the simplest case where A is diagonal, A\a) = a(a)\a), 
we find that A(a, S n ) = a(a) and 



p=0 



w 



(11) 



Next we treat the case where A = H m , where H m is 
one of the terms in the Hamiltonian, then 

i 00 on n 

(A) = (H m ) = ^EEE^H^n^i")- ( 12 ) 

a n=0 S n ' i=l 

and for each index sequence S n for the state |a), there is 
a sequence Sn+i, with H m as the last element. Defining 



(13) 



and considering Eq. (fiXty, the expectation value is ob- 
tained by counting the number N(m) of H m operators 
in the sequence, 



(H m ) 



(14) 



The energy is the negative sum of all operators H n 
and it therefore follows that 



E 



(15) 



and we see that the energy is proportional to the average 
length of the operator sequence. One may argue that 
this is a strange formula; it seems that one should be 
able to decrease the average string length by adding a 
positive constant to the Hamiltonian (thereby increasing 
the energy). One should even be able to set the energy 
and average string length to zero! The solution to this 
paradox lies in the infamous sign problem. Adding a 
(sufficiently large) positive constant to the Hamiltonian 
will make it impossible to keep the weight function pos- 
itive, and hence Eq. ([l5|) is no longer valid. In fact, for 
most models, a negative constant has to be added to the 
Hamiltonian in order to make W(a, S n ) positive definite. 
Eq. ( [l5|) then gives the energy including these constants. 

Other expectation values, such as products of 
imaginary-time dependent operators and static suscepti- 
bilities can also be easily evaluated with the SSE method. 
We refer to previous work for the expressions for these 
more complicated expectation values£3 

In order to make the simulation code efficient one as- 
sumes that only operator strings with a length shorter 
than L contribute to the trace. This is not necessary, 
but by automatically adjusting L so that the simulation 
will not reach strings longer than L in any reasonable 
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simulation time, only an exponentially small, completely 
undetectable, error is made. With the length of the op- 
erator string limited to L, (n — L) identity operators can 
introduced in an operator string of length n < L to make 
the string length fixed. For every original operator string 
of length n there are then (^) strings of length L, cor- 
responding to all possible insertions of the (L — n) unit 
operators. The inverse of this factor is included in the 
Taylor series expansion, where the summation now is over 
all index sequences of length L, denoted Sl- 

a S L ' i=l 

During the simulation the number of these extra unit 
operators will fluctuate, and hence the previous sum over 
n is implied. The advantage with fixing the length in this 
manner and introducing extra unit operators is that all 
moves can be defined in solely terms of exchanges, i.e. a 
set of one or several operators is exchanged for another 
set of the same number of operators. This simplifies the 
construction of an updating scheme that satisfies detailed 
balance. 



C. Application to the Heisenberg ferromagnet 

In this Section we give a detailed description of how the 
SSE method is applied to the Heisenberg ferromagnet. 
The Hamiltonian of this model including a magnetic field 
is given by 

H=-J^S i .S s(i) -hY^S x (i), (17) 

i,S(i) i 

where the coupling constant J > 0, h denotes the mag- 
netic field strength and 5(i) denote the nearest neighbors 
of site i (we count each interacting spin pair only once) . 
For reasons that will be explained later, the magnetic 
field is chosen in the x direction. The rectangular lat- 
tice has n = l\ x I2 sites, where l\ and I2 are the linear 
dimensions of the rectangular lattice. Throughout this 
work we use periodic boundary conditions, and there are 
therefore m = n bonds if l\ — 1, corresponding to the 
one-dimensional lattice, and m = 2n bonds if Zi > 1, 
corresponding to the two-dimensional case. 

For the purpose of the SSE updating scheme we intro- 
duce the following six operators: 



#0,0 


= I 


•Hi, 6 


= h 


H2,i 


= Ii 


#3,6 


= 2(Si(b)Sj{b) + jh) 






#5,6 


= S ](b) S k(b) + S j(b) S k(b) 



where i denotes a lattice site, b a bond, i(b) and j(b) 
are the sites connected by bond b, and I is an identity 
operator. The unit operators h and 7i, labeled by a site 
or bond, are introduced to simplify the updating scheme. 
An algorithm without these could also be formulated. 
The Hamiltonian can then be written as 

J m 

b=l 

hsr^, rr . 3Jm hn , „. 

--J2( H 2,i+ H ±,i) + — + ^- (19) 
»=i 

Note that the operator H = I is not considered 
a term of the Hamiltonian; it is the unit operator em- 
ployed to augment the operator strings to the fixed length 
L, as discussed in the previous Section. The constant 
in the definition of the type-3 operator has been intro- 
duced in order to make all its matrix elements positive (or 
zero), thereby making the weight function positive defi- 
nite. The introduced constants only shift the energy (and 
therefore also the mean length of the operator string). 

According to an ergodic updating scheme, the six dif- 
ferent kinds of operators in the operator sequence are 
interchanged in such a way that any contributing oper- 
ator string and basis state can be generated by a series 
of updates. Before going into the details of these pro- 
cedures, we will specify our basis states and the weight 
function, and introduce some notation. 

We will work in the S z basis \S* ,5f , . • . , S^), and the 
non-branching property is satisfied since 



flsl U) 


= H 3 \ IT) 


flsl TT) 


= 1 TT) 


H 3 \ li) 


= 1 II) 


Hi\ T) 


= 11) 


Hi\ i) 


= IT> 


h b \ TT) 


= H 5 \ ||) 


h b \ IT) 


= IU) 


H 5 \ Tl) 


= 1 IT). 



When non-zero, i.e., if \a(L)) — \a(0)), the weight 
function is 

W{a,S L )=(r{L-n)\r-\ f-J , (21) 

where is the total number of _fffc,i (i = 1, . . . , m) oper- 
ators in the operator string. The weight is always positive 
because all the matrix elements of the operators defined 
in Eq. (|2l]) are positive or zero. For convenience we now 
use a two- index notation for the operator-index sequence 
Sl', 

& = (*') ( fa ) ...( fa ) , (22) 
VJi/iVJ2/ 2 \3lJ L 
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|oc(0)> 

|a(5)> 
|a(4)> 

|a(3)> 

|a(2)> 

|a(l)> 

|a(0)> 



(2 4) 
(4 2) 
(4 1) 
(3 2) 
(0 0) 
(5 1) 



FIG. 2. Typical SSE configuration for a four-site system. 
The propagated states are shown on the left and the operator 
string to the right. The first operator- index denotes the type 
of operator (1 — 6) and the second is the bond on which it 
operates (1 — m), with m = 4 in this case. 



where k indicates the kind of operator (k G {0, • • • , 5}) 
and j indicates what site or bond the operator acts on. 

Next we introduce a random walk that satisfies the 
detailed-balance principle and covers the space {|a) (g> 
{S n , n = 0, . . . , L}}. A simulation is started from a ran- 
dom initial state |cc(0)) and an initial sequence of oper- 
ators (q) (q) • • • (g) l . We define six fundamental moves 
that together ensure that we cover the full configuration 
space as long as the magnetic field h is nonzero. In the 
absence of a magnetic field, additional moves have to be 
made, and these will be described later. The moves have 
to be carefully defined so that the detailed-balance prin- 
ciple is satisfied, generation of zero- weight configurations 
is not attempted, and so that they can be rapidly exe- 
cuted. 

To visualize the seemingly abstract SSE space, a typ- 
ical configuration for a one-dimensional four-site system 
is shown in Fig. ||. In the figure we can also see the close 
resemblance between the SSE method and the standard 
Euclidean path integralpJprmulations. This relation has 
been explored in detail.E2l Note that only the first state 
| £*(())) and the operator string have to be stored in mem- 
ory. All propagated states can be generated sequentially 
as needed, and the memory requirements for the method 
are therefore quite modest. 

The first kind of move changes the number of nonzero 
operators in the string by introducing the unit operator 
of type 1. Going through the operator string from the 
first to the last operator, attempts are made to exchange 
every (°) operator for a ( 6 ) operator (with b chosen 

randomly from {!,..., m}) and every (,) operator for a 



( Q ) operator. In an accepted move n and n\ in Eq. ( |2l| ) 
changes by ±1. The detailed-balance principle then re- 
quires that the probabilities P[Q <->(*)] of carrying 
out the replacements in the two different directions sat- 
isfy 



P{Q P ^(D P ] _ W{a,Sn +1 )p[Q p ^O p ] 



(23) 



where m again is the number of bonds and p[(g) *-* 
(j) ] denotes the a priori probability of the move being 
carried out, i.e., the probability before any accept/reject 
probability has been assigned. If the bond b is chosen 
at random for a move in the — > direction, the a priori 
probabilities satisfy p[Q p Q p ] = mp[(° Q ) p -> Q p ], 
since there are m different 1-operators, but only one 0- 
operator, i.e., every accepted transition from any of the b 
I-operators leads to the same 0-operator, but a transition 
from the 0-operator only leads to any I-operator for a 
given bond with probability 1/m. One can easily verify 
that we satisfy the detailed balance condition by choosing 



P 



P 



mf3J 
2{L - n) 

2(L-n+ 1) 
Jm/3 



(24) 



where a number greater than 1 on the right hand side 
should be interpreted as probability one. 

The second kind of move is very similar, but it ex- 
changes type-0 and type-2 operators. Again the string 
is sequentially searched for these two kinds of operators, 
and detailed balance is satisfied with the following ex- 
change probabilities: 



P 



P 



m(3h 
(L-n) 

(L-n + 1) 
hm[3 



(25) 



The third kind of move attempts to exchange type-1 
and type-3 operators. We can see that the weight func- 
tion is unaffected by this kind of move, and detailed bal- 
ance is satisfied, but there is a restriction. An attempt to 
exchange ( 6 ) for will result in a zero-weight func- 
tion if the spins on the two sites that are connected to 
bond b point in opposite directions. Hence one needs 
to check the spin configuration before attempting to ex- 
change a 1-operator for a 3-operator, while a 3-operator 
can always be exchanged for a I-operator. 

The fourth kind of move is slightly more complicated 
and involves exchanging pairs of 2- and 4-operators. The 
reason why we have to exchange pairs of operators is 
that the state of the system has to return to its original 
state when it has been propagated by the whole operator 
sequence; |a(0)) = \a(L)), and if a spin is flipped in an 
intermediate state it has to be flipped back at a later 
time. Hence we attempt to make exchanges of the form 
Q)pQ)q (T) P (T)q- But we nave to be careful, since this 
move results in the spin at site i being flipped in all states 
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\a(p)), . . . , \a(q — 1)}. In case p > q, then |a(0)) will be 
affected, and we have to flip the spin at site i in the state 
|a(0)) that we store in the memory. Also, if there is an 
operator of type (?) or y) with p < r < q, then this 
exchange would lead to a zero weight function and the 
move must not be made. In order to make an efficient up- 
dating scheme the operator sequence is divided up into n 
subsequences, each containing the necessary information 
to make the above exchange for site i <E {1, . . . n}. 

It is easiest to actually store three lists: Ti,Pi and 
Fi, that contain information about the ith site. The list 
Ti = {ti, . . . , iL(i)}> with t £ {2, 4} contains all 2- and 4- 
operators from the full operator string operating on site 
i. The list Pi = {p%, ■ ■ ■ ,PlU)}i contains the positions p 
of these operators in the full operator string. The final 
list is Ft — {/i, . . . , where fj £ {0, 1} indicates if 

there are any operators of kind 3 or 5 between the po- 
sitions pj and Pj+i- The move then simply consists of 
choosing a tj with the corresponding fj = 0. Thereafter 
one moves forward in the F list until the first nonzero fi is 
encountered. Then one of the coordinates k (j < k < I) is 
chosen and if tj = tk the move is accepted, but if tj =/= t% 
the move is rejected, since the two operators then are of 
a different kind. If tj ^ tk one can, however, permute 
the two operators. Again the weight function in unaf- 
fected by this move, and the detailed-balance principle is 
automatically satisfied. 

The fifth move involves another pair exchange of the 

form (l)p(l) q ~ (i)p(iV A S ain we have t0 be careful 
since this move flips the two spins connected to bond b 
in all states \a(p)), . . . , \a(q — I)}, and the weight func- 
tion will vanish if there is an operator of type ( fe ,) or 

( b ,) , with p < r < q, where the bond b' is connected 
to either of the sites that bond b is connected to. The 
full operator sequence is therefore again divided up into 
substrings containing the operators that act on a par- 
ticular bond and the sites that it is connected to. One 
can make this kind of move in two different substrings 
independently of each other, as long as the two bonds do 
not connect, i.e. as long as they do not have any site in 
common. The n bonds on a one-dimensional lattice can 
be partitioned in all the bonds that start at an odd site, 
or all the bonds that start on an even site. Hence the full 
operator string is split up two times, each time into n/2 
sub-strings, which can be updated independently. A two- 
dimensional lattice has to be divided into four separate 
partitions; see Fig. ||. 

The updating of the n/2 sublists is very similar to 
the updating scheme for the fourth move. Four lists 
are stored this time. The list Tj = {ti, . . . , t^a)}, 
with t G {1>5} contains all 1- and 5-operators from 
the full operator string operating on bond b. The list 
Pb = {pit ■ ■ • jPlU)} again contains the positions p of 
these operators in the full operator string. The list 
^6 = {/i> • • • j /L(i)}j where fj 6 0,1 indicates if there 
are any operators of kind 3 or 5 operating on a bond 
connected to bond b between the positions pj and Pj+\- 
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FIG. 3. The four different bond partitions of a 6 x 6 lattice 
with periodic boundary conditions. 



The final list is Sb = {s±, . . . , s^m}, where Sj = 1 in- 
dicates that the spins connected to bond b are aligned 
in state \a(pi)), and Si = indicates that they are anti- 
aligned. The move is identical to the procedure described 
for the fourth move, except that now the move is canceled 
if tj =/= tk or Sj = 1. 

The sixth move is the most complicated one, because 
it involves exchanges between three different kinds of op- 
erators. It is this move that makes this simulation par- 
ticularly efficient, since it can generate configurations of 
non-zero winding number. 

The following exchanges are attempted: (?) (?) «-» 

(b) P (h) 9 ' wnere H an d *2 are the two sites connected to 
bond b. Again the exchange involves spin flips, and the 
lattice is divided up into partitions as in move 5. Four 
lists are again used, where X& includes all l-,4- and 5- 
operators, Pb gives the position of the operators in the 
full sequence, Fb indicates whether an operator of type 3 
acts on the two sites and Sb keeps track of the spin con- 
figuration. In the same fashion as above two operators, 
tj and tk, are chosen, and if {tj, tk} = {1 , 5} an attempt 
to change the two operators to tj — tk — 4 is made, if 
the spin configuration allows for this. If tj = tk = 4 an 
attempt to change the two operators to {tj, tk} = {1,5} 
is made. The detailed-balance principle is satisfied if 



1 

b) \bj 

p \ / q 



*l/p V2, q 



L s p \ q 



1\ /5 



b \b, 

' p \ / q 



= 2h 2 /J 2 (26) 
= J 2 /2h 2 . (27) 



If the last move is excluded one can see that all possible 
operator sequences are not sampled if periodic bound- 
ary conditions are used. To visualize this, we consider a 
four-site system for which is is easy to realize that the 
simple string iH^^Flb, 3H5, 4 cannot be reached. This 
is an example of a configuration with a non-zero winding 
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number (provided that the operator string operates on 
an appropriate state). In order to make this configura- 
tion possible we need to introduce a "ring" move. If we 
picture each 5-operator as connecting the two sites it acts 
on, we can see that the configuration H§_ 1^5,2-^5, 3^5, 4 
creates a ring around the system. The ring move is ac- 
complished as follows: starting at a random point in the 
operator sequence for a n-site system a search for a set 
of n/2 different 5-operators is made. If such a set is 
found, an attempt to exchange it with its complementary 
set is made. An example for the 4-site system would be 
#5, life 3 — > H5 2H5. 4- If the move is successful the wind- 
ing number has changed. Whether or not the move can 
be carried out of course depends on constraints imposed 
by the states. As the system size is increased beyond 
about L = 16 it becomes virtually impossible to perform 
ring moves. 

In a one-dimensional system the ring move is the only 
move that has to be added if the sixth move is not per- 
formed. In two dimensions we can picture a ring around 
a small part of the system. We can, for example, draw a 
ring around a 4 x 4 square and we realize that we can- 
not reach a configuration consisting of the 14 5-operators 
that connect these sites. We therefore introduce a pla- 
quette move. This move changes the 5-operators that act 
on a single plaquette, which for a two-dimensional square 
lattice is the smallest square in the lattice. A plaquette 
move that involves arbitrarily large parts of the system, 
for example the 4x4 square considered above, can be 
reached with these fundamental plaquette moves. The 
move is identical to the one-dimensional ring move for a 
4-site ring. 

In higher dimensions, when using periodic boundary 
conditions, we also have to perform a direct general- 
ization of the one-dimensional ring move. In a two- 
dimensional system with periodic boundary conditions 
we can picture making rings around the system in both 
spatial directions, and in three dimensions we could make 
the move in all three spatial directions. Apart from per- 
forming the move in all spatial directions it is identical 
to the ID ring move, and for a linear system size larger 
than about L = 16 ring moves are no longer accepted. 
Note that this ring move can not be accomplished by the 
above plaquette moves. 

Without the external field in the x-direction, the 
Hamiltonian conserves the total magnetization in the z- 
direction; M z = Sf, and the operator string updates 
alone therefore sample only within a fixed magnetization 
sector. In order to sample in the grand canonical ensem- 
ble, one then has to carry out an update that changes the 
magnetization of the state \a) in Eq. (Jfy. A spin Sf in \a) 
can be flipped, without change in the weight, if there are 
no operator L) or (jj) in Sl with b a bond connected to 
spin i. The probability of this being the case approaches 
zero rapidly as T is lowered, and hence this update can 
be carried out only at relatively high temperatures. With 
an external field in the x direction the magnetization is 
no longer conserved and the simulation is automatically 



grand canonical. It is nevertheless useful to carry out the 
single-spin flips at high temperatures. 

With the field present in two dimensions, the sixth 
move also makes plaquette and ring moves unnecessary. 
This can be understood since the sixth move introduces 
single 5-operators, and through a series of moves any 
configuration of 5-operators can be generated. Had we 
chosen the field in the z-direction this would not have 
been possible. Having the field in the x-direction there- 
fore makes the simulation particularly efficient by intro- 
ducing single-spin flipping operators. Strictly speaking 
we do not need to make the fifth move either, but espe- 
cially at small fields it may be good to include this move 
to shorten the auto-correlation time and speed up the 
thermalization of the system. 

Another advantage of choosing the n-direction for the 
field is that we can then easily measure spin-spin corre- 
lations both parallel and perpendicular to the field. The 
perpendicular correlation functions are diagonal in our 
chosen basis, and can be evaluated according to Eq. (O). 
The parallel correlations involve expectation values of 
products of the number of field operators in Sl as dis- 
cussed in Refs. [l?]. We are particularly interested in the 
perpendicular correlations for this model, since they are 
the ones required to calculate the relaxation rate 1/T\. 
The magnetization M = (Sf + S^)/2 is according to 
Eq. (|l4|) proportional to the average total number of field 
operators H± ti in the sequence; M = (n 4 ) / '(N f3h). 

A complete Monte Carlo step for the two-dimensional 
lattice, as used in this work, consists of 

1. Move 1, Qj) <-> (l) p , attempted at all positions 
p = 1, . . .L in the sequence Sl- 

2. Move 2, <-» (J , attempted at all positions 
p = X, ... L in Sl- 

3. Move 3, Q) <-» (j) j attempted at all positions 
p = X, ... L in Sl- 

4. The sequence is split up into n sub-sequences and 
in each sequence move 4,(?) (^) <-> (J , (J is 

attempted, whereafter the full operator sequence is 
restored. 

5. The operator string is split up 4 times into n/2 
different sub-sequences. In each sub-sequence the 
fifth move, (l) p (l) q «-► (l) p (l) q , is then attempted 
a number of times of the order of the length of 
the subsequence. The updated sub-sequences are 
recombined into the full string. 

6. The operator string is split up 4 times into n/2 
different sub-sequences and in each sub-sequence 




, is at- 



tempted. The sub-sequences are recombined to the 
full string. 
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7. All spins that are not acted on by any interaction 
operators (.) or (,) are flipped with probability 
1/2. 
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FIG. 4. Autocorrelation function for the x-component of 
the magnetization of 8 x 8 lattices at three different strengths 
of the external field. The upper panel is for temperature 
T/J = 1, and the lower panel for T/J = 1/4. 
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FIG. 5. The time dependence of the number of field opera- 
tors in the operator string at T/J — 1 and two different field 
strengths. 

In principle it is possible to combine moves 1-3 into a 
single move with several "branches" . For simplicity we 
have not discussed this (only marginally) more compli- 
cated approach here. 

A full simulation consists of a number of equilibration 
steps, during which the maximum length of the string is 
automatically increased as the string grows. The length 
of the string very quickly reaches its equilibrium, and af- 
ter the equilibration steps are done the maximum length 
is fixed (as described in the previous section) and mea- 
surements are carried out at even intervals. 

The present updating scheme is very efficient provided 
that the density of single-spin flipping operators Ha in 



the sequence is not too small. In practice, we have found 
that simulations for h/J > 0.05 deliver accurate results 
without too much effort, independently of T. In Figure 
H we show some results for the auto-correlation function 
of the magnetization, defined according to 



Cm At) 



{M x (i)M x {i + i))-{M x f 
(MJ) - (M s >2 



(28) 



where M x (i) is the value of the magnetization estimator, 
714/ \N(3h), at the ith Monte Carlo step. We show results 
for 8 x 8 lattices at two different temperatures and three 
different field strengths. Note that the long-time cor- 
relations decrease with decreasing temperature. This is 
related to the fact that the distribution of magnetization 
values becomes broader. For example, for h/J = 0.05 
at T/J = 1, the magnetization essentially fluctuates be- 
tween two values, corresponding to configurations with 
or 2 field operators, as shown in Figure ||. The sys- 
tem spends long times in states with no field operators 
and occasionally short times with 2 operators. These two 
time scales are reflected in the autocorrelation function, 
which shows a rapid decay at short times, but a very 
slow decay at longer times (the asymptotic exponential 
decay time is several hundred Monte Carlo steps). At 
the lower temperature the short-time behavior shows a 
less rapid decay, but also a faster asymptotic decay, as 
the fluctuations in the distribution of the number of field 
operators TI4 is now broader, and the likelihood of a fluc- 
tuation by ±2 in the simulation is higher. At h/J = 0.2, 
the time dependence of 714 shows fluctuations on much 
shorter time-scales, and the autocorrelation function ac- 
cordingly decays considerably faster. 



III. THE MAXIMUM ENTROPY TECHNIQUE 

It is notoriously difficult to obtain dynamic informa- 
tion from quantum Monte Carlo simulations. One of the 
presently most common and-oromising methods is the 
maximum entropy technique,Ej which we used to calcu- 
late the spin-lattice relaxation rate. In this section we 
will discuss various ways to use the method, compare 
results to exact diagonalization and show some of the 
limitations of the procedures. 

The nuclear magnetic resonance spin-lattice relaxation 
rate is given by 123 



2 S(q,uj N ), 



(29) 



where A q is the Fourier transform of the hypcrfme cou- 
pling. The resonance frequency wjy is so small compared 
to other energy scales that we consider the limit ui — > 0. 
From now on we assume that A q = 1. The dynamic 
structure factor S(q,u), which measures transverse spin 
correlations, can be obtained from the imaginary time 
dependent correlation function 



8 



C(q, r)=Yl T ) = E ^ r (sf(r)S! +r {0)) (30) 

r r 

by inverting the expression 
£C(?,r)= [ <L;Y / S(q,cj)e-^ = C(r = 0,T). (31) 

Calculating the local imaginary time correlation function 
is therefore in principle sufficient for determining the re- 
laxation rate, after performing a maximum entropy in- 
version. For reasons that will become clear below, we 
did, however, also measure the spatial dependence of the 
correlation function. 

Using conservation of momentum it is possible to ob- 
tain finite temperature exact diagonalization results for 
system sizes up to 4 x 4. Below we will show diagonaliza- 
tion data compared to Monte Carlo + maximum entropy 
results. 

First a few words concerning the default model needed 
in the maximum entropy method. All data shown in this 
paper have been generated using a flat default model as 
defining zero entropy. Much has been written about how 
to choose a good default modelpd but from our expe- 
rience it seems that unless something very specific and 
very accurate is known about the solution (in which case 
it probably is unnecessary to use Monte Carlo + max- 
imum entropy), one should use a flat default model. If 
only general features, obtained from, for example, per- 
turbation theory or some analytical mean field solution, 
are known about the function, then one introduces a bias 
towards one approximate model by placing peaks at fre- 
quencies that are only approximate. If the answer is very 
dependent on what default model one uses, then the accu- 
racy of the results should be viewed with caution. Using 
good data and a flat default model is optimal in the sense 
that no prior bias has been used to construct features in 
the spectral function; the default model only has a reg- 
ularizing effect. We worked extensively with Gaussian 
default models, where we determined all the parameters 
of the default model through sum rules that could be 
calculated in the Monte Carlo simulation. In some cases 
this worked fairly well, but in those cases a flat default 
model worked almost equally well. In other cases the 
maximum entropy solution would be very close to the 
Gaussian default model, but far from the exact diagonal- 
ization results. This clearly shows that the problem is 
ill-posed. In these instances it was our experience that 
a flat default model resulted in better agreement with 
exact diagonalization. 

We now turn to actual comparisons with exact diago- 
nalization results. At low temperatures the spectral func- 
tion is dominated by a delta function for each magnon, 
and this extreme limit will first be considered to illustrate 
some weaknesses and strengths of the maximum entropy 
solution. In Fig. ^ we show results for a 4 x 4 system 
at inverse temperature T/J = 0.25 and field strength 
h/J = 0.25. As usual, the exact diagonalization result 
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FIG. 6. Maximum entropy compared to exact diagonaliza- 
tion for a 4 x 4 system at T/J = 0.25 and h/J = 0.25. The 
exact results are represented by the solid-line histogram. The 
maximum entropy method is applied to each momentum sep- 
arately (dotted curve) , to the average of all momenta (dashed 
curve) and to the average of all momenta, except the q — (0, 0) 
(long-dashed curve). 

really only consists of a series of delta peaks, but the spec- 
trum has been divided into bins of finite width in order to 
illustrate the result more clearly. The first peak is at ex- 
actly io — h/J = 0.25 and is the response of the q = (0, 0) 
momentum. This peak will remain a delta function even 
at finite temperatures since the spin-correlation function 
then will contain the total spin raising and lowering oper- 
ators, and hence transitions can occur only between levels 
within the same spin multiplet. These levels are all sep- 
arated by the Zeeman splitting, which is independent of 
temperature. Therefore the delta peak at the Zeeman 
energy will not get broadened by temperature. For finite 
q-values there will be transitions between different spin 
multiplcts causing transition energies that differ from the 
Zeeman energy. 

In order to be able to directly study the different mo- 
menta we have measured the full spatial dependence of 
the correlation function and we can therefore either ana- 
lytically continue each momentum and then average the 
spectrum, or first average the correlation function and 
then do the analytic continuation. In Fig. [6] both meth- 
ods are demonstrated. The dotted curve shows the result 
when analytically continuing each momentum separately. 
We notice two important features. First, the maximum 
entropy method resolves each peak, and secondly, the 
resolution of the peaks decreases with frequency, as can 
be expected, since the factor e~™ makes the inversion 
insensitive to high frequency features. The dashed curve 
shows the result when the average over all momenta is 
continued. We notice immediately that the maximum 
entropy method has difficulty in resolving more than one 
peak and tends to smear the result out. Focusing on the 
limit u^Owe notice that because of the q = (0, 0) peak 
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FIG. 7. Maximum entropy results compared to exact di- 
agonalization data for a 4 x 4 system at T/J = 2.0 and 
h/J = 0.25. The exact results are represented by the 
solid-line histogram. The maximum entropy method is 
applied to each momentum separately but not including 
q = (0,0) (dotted curve), to each momentum separately in- 
cluding (0,0) (dashed curve), to the average of all momenta 
not including (0, 0) (long-dashed curve) and to the average of 
all momenta including (0,0) (dot-dashed curve). 

at u> = h/J = 0.25, an estimate of 1/T\ would be too 
high, since the weight of the q — (0, 0) peak is smeared 
out all the way to ui = 0. To remedy this effect we can 
average over all momenta, except q — (0, 0). The results 
are shown in the dot-dashed curve, and we see that this 
certainly affects the solution close to w = 0, where there 
is now correctly, no weight. 

We are, however, interested in results for large system 
sizes, for which the spectral function should be smooth 
on a reasonably fine frequency scale. We cannot com- 
pare our results to diagonalization results for more than 
4x4 sites. We can, however, study this system at higher 
temperatures, for which the spectral function is fairly 
smooth. In Fig. fj] the same parameter values as in the 
previous figure are shown, but at a higher temperature, 
T/J = 2.0. We see that the diagonalization results are 
much smoother and we can analyze the maximum en- 
tropy results. The dashed curve shows all momenta sep- 
arately continued. l/7\ is grossly overestimated since 
the q = (0, 0) peak is smeared out at this higher tem- 
perature. Removing the q = (0, 0) peak gives a much 
better estimate; see the dotted curve. Continuing the 
average of all momenta again places too much weight at 
low frequencies, because of the q = (0, 0) peak, as can be 
seen from the dot dashed curve. Removing the q = (0, 0) 
momentum from the average improves the result a great 
deal; see the long dashed curve. 

We found that the optimal use of the maximum en- 
tropy method in our case was to continue the average of 
all momenta, excluding the q = (0, 0) momentum. For 
this momentum we know the exact analytic result, and 



including it in the continuation leads to an over-estimate 
of \/T\. The reason for not continuing each separate 
momentum, which worked miraculously well in Fig. ^j, is 
that at intermediate temperatures this method also leads 
to an overestimation of 1 /T\ , since peaks close to the ori- 
gin will be smeared out, and much of their weight will be 
incorrectly placed at ui — 0. 

To conclude this Section, we have found that the max- 
imum entropy method is a useful method to obtain real 
time data from imaginary time quantum Monte Carlo 
data for the Heisenberg ferromagnet. We believe that 
the specific recipe for how to use the method probably 
varies from system to system, but for the case we have 
considered here we found the most useful default model 
to be flat, and we also found that we had to separate out 
the q = (0, 0) momentum before carrying out the analytic 
continuation. 

Similar calculations of NMR relaxation rates have pre- 
viously beep-jcarried out also for various quantum antifer- 
romagnets.Ea In that case the differences between maxi- 
mum entropy calculations based on momentum and real- 
space correlations were less pronounced than what we 
have found here. The ferromagnetic case appears to be 
more difficult since the field induces additional structure 
in the frequency dependence. Also, at low temperatures 
independent magnons are exact excitations of the ferro- 
magnet, which causes further sharp features in the spec- 
tral function. 



IV. RESULTS 
A. Magnetization 

The Heisenberg model is one of the basic non-trivial 
models of a quantum ferromagnet. A detailed knowledge 
of the field dependence of the magnetization is there- 
fore of interest for comparisons both with analytical and 
experimental results. In this Section we will show the 
general field dependence of the magnetization for the 
two-dimensional model. We will also compare our Monte 
Carlo results to exact diagonalization results and show 
the nature of the finite size corrections. We also compare 
the results to recent experimental results for the v = 1 
quantum Hall state. 

In Fig. ^| we show a comparison of exact diagonaliza- 
tion and Monte Carlo data for a 4 x 4 system. We have 
verified the agreement to a relative error of 10 -4 . To this 
degree of accuracy the results for the largest system sizes 
shown (16 x 16 and 32 x 32) have also converged. We did 
additional test runs for systems of size 64 x 64, and the 
results were within error bars of the above results. From 
Fig. H we notice that the finite size effects are largest 
at intermediate temperatures, which is easily understood 
since in the limit of very low temperatures the magneti- 
zation will be fully saturated independently of the lattice 
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size, while in the high temperature limit the magnetiza- 
tion will vanish independently of the system size. 
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FIG. 8. Monte Carlo results for the magnetization vs tem- 
perature for LxL systems with L = 4, 8, 16 and 32, and exact 
diagonalization results for L = 4. The magnetic field h — 0.1. 




FIG. 9. Magnetization vs temperature for the 2D Heisen- 
berg ferromagnet in magnetic fields of strengths h=1.00, 0.80, 
0.60, 0.40, 0.20,0.10 and 0.05, from top to bottom. The re- 
sults were calculated using lattices sufficiently large (typically 
32 x 32) to eliminate finite-size effects almost completely. 

In Fig. |o| we present a plot of the magnetization as 
a function of temperature for a range of field strengths. 
We have tried to scale the data as a function of T/h a for 
different values of the exponent a. Such scaling does not 
collapse the data onto a single curve, however. Analytical 
calculations have also suggested that there is no scaling 
in T and h for this model.0 

We also show a comparison of recent magnetoabsorp- 
tion measurementsn performed on the v = 1 quantum 
Hall state compared to quantum Monte Carlo data and 
Schwinger boson calculations in Fig. |l(J There are no 
adjustable parameters in this comparison, since J can be 
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FIG. 10. Magnetization curve for the v — 1 quantum 
Hall ferromagnet. Comparison of Schwinger boson [SU(iV) 
and O(N) symmetric versions of the theory], quantum Monte 
Carlo and experimental data (Ref. 2). The calculations were 
carried out at field strength h/J = 0.32. 

calculated exactly for the zero- width welH using the ex- 
actly known spin wave dispersion for the quantum Hall 
ferromagnet. and estimated for the actual experimen- 
tal systemci The excellent agreement shows that this 
itinerant ferromagnet is well described by an effective. 
Heisenberg model. As discussed in our previous papem 
the comparison with the analytical Schwinger boson so- 
lutions showed that the 1/N corrections to the O(A) 
model agrees well with Monte Carlo and experimental re- 
sults at moderate and intermediate temperatures, while 
the low-temperature result was better reproduced by the 
SU(A) model. For a more comprehensive discussion of 
the Schwinger boson results and the experimental data, 
we refer to our previous paper 



B. Spin-lattice relaxation rate 



In Fig. 
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we compare mean-field resultsS'i for the 
nuclear magnetic relaxation rate 1/T±, to numerical re- 
sults obtained by analytic continuation of imaginary 
time Monte Carlo data, as described previously. As we 
have discussed, the maximum entropy results have to be 
viewed with some caution. The error bars are obtained 
using the bootstrap technique.o They do not strictly cor- 
respond to a statistical likelihood that the estimated re- 
laxation rate is within error bars of the true relaxation 
rate, since there is also an unknown systematic error due 
to the bias of the maximum entropy procedure. The er- 
ror bars do contain information about how sensitive the 
results are to the variations in the MC data. 

We notice a fairly good agreement between the ana- 
lytic and numerical results and, as in the case for the 
magnetization^ the O(A) theory appears to be some- 
what closer to the numerical results. Notice that this 
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FIG. 11. Monte Carlo and maximum entropy numerical 
results (circles) and Schwinger boson results [SU(N) mean 
field (dashed curve) and O(N) mean field (solid curve)] for 
the relaxation rate 1/Ti for h/J = 0.10. 



maximum entropy method, and compared also these re- 
sults with the Schwinger boson results. The role of the 
maximum entropy default model was discussed, and for 
the 2D Heisenberg model we argued that a flat default 
model works best. We further discussed advantages and 
problems arising when applying the maximum entropy 
procedure to imaginary time correlation functions in real 
space and momentum space. 
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time the numerical results do not lie between the O(N) 
and SU(iV) solutions, which was the case for the magne- 
tization. At zero temperature the relaxation rate is zero, 
since no spins flip in the ordered ferromagnet and there- 
fore the nuclear spins cannot relax. At low temperatures 
the rate is activated and caused by thermally activated 
spin waves. 

Not enough experimental data is available for the re- 
laxation rate to make a comparison with the above nu- 
merical results. Because results for the Heisenberg model 
agreed very well with experimentally measured magneti- 
zation, it would be of interest to compare the relaxation 
rate, to see if the Heisenberg model also captures the 
correct dynamic behavior of the v = 1 quantum Hall 
ferromagnet. 



V. SUMMARY AND CONCLUSION 

We have described an approximation-free quantum 
Monte Carlo technique and applied it to a two- 
dimensional ferromagnet in a magnetic field. We have 
shown that applying the field in the transverse direc- 
tion causes the simulation to become free of system- 
atic errors although only local Monte Carlo moves are 
made. The transverse field also enables easy access to 
transverse spin correlation functions. We have calcu- 
lated the temperature dependence of the magnetization 
for a range of field strengths and discussed finite-size ef- 
fects. Some results have previously been compared to 
measurements of the v — 1 quantum Hall state and an- 
alytical Schwinger Boson calculations.!] The high rela- 
tive accuracy (10 -4 ) of the Monte Carlo results made it 
possible to make statements about the relative merits of 
the SU(iV) and O(iV) solutions. We have, in addition, 
calculated the spin-lattice relaxation rate 1/Ti using a 
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